North Sea SMS Keyrun Supplemental Plots
Valerio Bartolino, Sarah Gaichas
10/14/2020
# load the species list (checked that they are the same for the 2020 and 2017 runs)
spp <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/NorthSeaKeyRun_2017/species_names.in", nrows=27)
spp <- gsub("_","",spp$V1)
spp <- data.frame(Species.n=1:length(spp), stkName=spp)Observed Data vs SMS Predicted
Catch
2020 SMS KeyRun
# plot Obs vs Pred catch (here only 2020 run)
dat <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/KeyRunComparisons/NorthSeaSMS2020/summary_table_raw.out", sep="", header=T , na.strings="", stringsAsFactors=F)
#dat[1:3,]
dat <- left_join(dat,spp) %>%
mutate(Yield = as.numeric(as.character(Yield)),
Yield.hat = as.numeric(as.character(Yield.hat)),
Year = as.numeric(as.character(Year))) %>%
gather("Source","Yield",8:9)
#postscript(paste(dirFigs,"catch_obsVSpre.ps",sep="/"))
ggplot(dat) +
geom_point(aes(Year,Yield,col=Source)) +
geom_line(aes(Year,Yield,col=Source,lty=Source)) +
facet_wrap(~stkName, scales="free_y")+
xlim(1973,2020)+
theme_tufte() +
theme(legend.position="bottom")
#dev.off()2017 SMS KeyRun
# plot Obs vs Pred catch (here only 2020 run)
dat <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/NorthSeaKeyRun_2017/summary_table_raw.out", sep="", header=T , na.strings="", stringsAsFactors=F)
#dat[1:3,]
dat <- left_join(dat,spp) %>%
mutate(Yield = as.numeric(as.character(Yield)),
Yield.hat = as.numeric(as.character(Yield.hat)),
Year = as.numeric(as.character(Year))) %>%
gather("Source","Yield",8:9)
#postscript(paste(dirFigs,"catch_obsVSpre.ps",sep="/"))
ggplot(dat) +
geom_point(aes(Year,Yield,col=Source)) +
geom_line(aes(Year,Yield,col=Source,lty=Source)) +
facet_wrap(~stkName, scales="free_y")+
xlim(1973,2020)+
theme_tufte() +
theme(legend.position="bottom")
#dev.off()Diet proportion 2020
# plot Obs vs Pred stomachs (here only 2020 run)
dat <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/KeyRunComparisons/NorthSeaSMS2020/summary_stom.out", header=T)
#dat[1:4,]
sppPred <- spp %>% mutate(Predator.no = Species.n) %>% rename(PredName=stkName)
sppPrey <- spp %>% mutate(Prey.no = Species.n) %>% rename(PreyName=stkName)
dat <- dat %>%
left_join(sppPred) %>%
mutate(Species.n = NULL) %>%
left_join(sppPrey) %>%
mutate(PreyName = as.character(PreyName),
PreyName = ifelse(is.na(PreyName), "otherfood", PreyName))
tmp <- dat %>%
group_by(Year,PredName,PreyName) %>%
summarise(stomcon = sum(stomcon,na.rm=T),
stomcon.hat = sum(stomcon.hat,na.rm=T))
# pie plot Obs and Pred
preys <- as.character(unique(tmp$PreyName))
preds <- as.character(unique(tmp$PredName))
#plList <- vector("list",length(preds))
nb.cols <- 14
colPalette <- colorRampPalette(brewer.pal(12, "Paired"))(nb.cols)
plist2 = lapply(split(tmp, tmp$PredName), function(d) {
#for(i in 1:length(preds)){
ggplot(d %>% gather("source","stomcon",4:5)) +
geom_bar(aes(x="",stomcon,fill=PreyName), width = 1, stat = "identity", position="fill") +
scale_fill_manual(values=colPalette) +
facet_grid(source~Year) +
#ggtitle(preds[[i]]) +
coord_polar("y")
})
# for(i in 1:length(plList)){
# postscript(paste(dirFigs,"/stomachs_fitting_",preds[i],"_v1.ps",sep=""))
# print(plList[[i]])
# dev.off()
# }
# same thing as barplot plot Obs and Pred
#for(i in 1:length(preds)){
plist3 = lapply(split(tmp, tmp$PredName), function(d) {
ggplot(d) +
geom_col(aes(Year-0.2,stomcon,fill=PreyName), width=0.35, stat="identity", position="stack") +
geom_col(aes(Year+0.2,stomcon.hat,fill=PreyName), width=0.35, stat="identity", position="stack") +
scale_fill_manual(values=colPalette) +
#ggtitle(preds[[i]]) +
xlab("Year") +
theme_bw()
})
# for(i in 1:length(plList)){
# postscript(paste(dirFigs,"/stomachs_fitting_",preds[i],"_v2.ps",sep=""))
# print(plList[[i]])
# dev.off()
# }Pie charts
for(i in 1:length(preds)) {
cat(" \n####", preds[i]," \n")
print(plist2[preds[i]])
cat(" \n")
}Cod
$Cod 
Fulmar
$Fulmar 
Gannet
$Gannet 
GBB.Gull
$GBB.Gull 
Guillemot
$Guillemot 
Haddock
$Haddock 
Her.Gull
$Her.Gull 
Kittiwake
$Kittiwake 
Mackerel
$Mackerel 
Puffin
$Puffin 
Razorbill
$Razorbill 
Saithe
$Saithe 
Whiting
$Whiting 
Greyseal
$Greyseal 
H.porpoise
$H.porpoise 
N.horsemac
$N.horsemac 
A.radiata
$A.radiata 
G.gurnards
$G.gurnards 
W.horsemac
$W.horsemac 
Hake
$Hake 
Bar charts
for(i in 1:length(preds)) {
cat(" \n####", preds[i]," \n")
print(plist3[preds[i]])
cat(" \n")
}Cod
$Cod 
Fulmar
$Fulmar 
Gannet
$Gannet 
GBB.Gull
$GBB.Gull 
Guillemot
$Guillemot 
Haddock
$Haddock 
Her.Gull
$Her.Gull 
Kittiwake
$Kittiwake 
Mackerel
$Mackerel 
Puffin
$Puffin 
Razorbill
$Razorbill 
Saithe
$Saithe 
Whiting
$Whiting 
Greyseal
$Greyseal 
H.porpoise
$H.porpoise 
N.horsemac
$N.horsemac 
A.radiata
$A.radiata 
G.gurnards
$G.gurnards 
W.horsemac
$W.horsemac 
Hake
$Hake 
SMS output comparisons
Numbers at age in Q1
# plot Num@age from 2017 and 2020 runs
dat17 <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/NorthSeaKeyRun_2017/summary.out", sep="", header=T , na.strings="", stringsAsFactors=F)
dat20 <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/KeyRunComparisons/NorthSeaSMS2020/summary.out", sep="", header=T , na.strings="", stringsAsFactors=F)
dat17 <- left_join(dat17,spp) %>% mutate(run="r17")
dat20 <- left_join(dat20,spp) %>% mutate(run="r20")
dat <- bind_rows(dat17,dat20)
# select Herring in Q1
#i <- spp$stkName[21] # 21 is herring
#postscript(paste(dirFigs,"/numAtAge_",i,".ps",sep=""))
plist = lapply(split(dat, dat$stkName), function(d) {
ggplot(d %>% filter(Quarter == 1)) +
geom_point(aes(Year,N, col=run)) +
geom_line(aes(Year,N, col=run)) +
facet_wrap(~Age, scale="free_y")+
theme_tufte() +
theme(legend.position="bottom")
})
#dev.off()Cod
plist$Cod
Haddock
plist$Haddock
Herring
plist$Herring
Mackerel
plist$Mackerel
N. Sandeel
plist$'N.sandeel'
Norway pout
plist$'Nor.pout'
Saithe
plist$'Saithe'
S. Sandeel
plist$'S.sandeel'
Sprat
plist$Sprat
Whiting
plist$Whiting
Modeled Q1 predator diet
2020 SMS KeyRun
# plot predators model diet (here only 2020 run)
dat20 <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/KeyRunComparisons/NorthSeaSMS2020/who_eats_whom_level2.csv", sep=",", header=T)
#dat[1:4,]
# aggregate birds
datQ20 <- dat20 %>%
mutate(Predator = as.character(Predator),
Predator = ifelse(Predator %in% c("Fulmar","Gannet","GBB. Gull","Gillemot","Her. Gull","Kittiwake","Puffin","Razorbill"), "Birds", Predator)) %>%
group_by(Year,Quarter,Predator,Prey) %>%
summarise(eatenW = sum(eatenW, na.rm=T))
# All predators diet in Q1
#postscript(paste(dirFigs,"diet_byPredator_Q1.ps",sep="/"))
ggplot(datQ20 %>% filter(Quarter == 1)) +
geom_bar(aes(Year,eatenW,fill=Prey), stat="identity", position="fill") +
facet_wrap(~Predator)+
xlim(1973,2020)+
theme(legend.position="bottom")
#dev.off()2017 SMS KeyRun
# plot predators model diet **check to make sure this is 2017 output**
dat17 <- read.table("https://raw.githubusercontent.com/ices-eg/wg_WGSAM/master/Input_Output/Output/WhoEatsWhom/who_eats_whom_level2.csv", sep=",", header=T)
#dat[1:4,]
# aggregate birds
datQ17 <- dat17 %>%
mutate(Predator = as.character(Predator),
Predator = ifelse(Predator %in% c("Fulmar","Gannet","GBB. Gull","Gillemot","Her. Gull","Kittiwake","Puffin","Razorbill"), "Birds", Predator)) %>%
group_by(Year,Quarter,Predator,Prey) %>%
summarise(eatenW = sum(eatenW, na.rm=T))
# All predators diet in Q1
#postscript(paste(dirFigs,"diet_byPredator_Q1.ps",sep="/"))
ggplot(datQ17 %>% filter(Quarter == 1)) +
geom_bar(aes(Year,eatenW,fill=Prey), stat="identity", position="fill") +
facet_wrap(~Predator)+
xlim(1973,2020)+
theme(legend.position="bottom")
#dev.off()Modeled predators contributing to predation by prey
2020 SMS KeyRun
# ------------------------
# plot predators contributing to predation by prey
preys <- as.character(unique(datQ20$Prey))
plList <- vector("list",length(preys))
nb.cols <- 14
colPalette <- colorRampPalette(brewer.pal(12, "Paired"))(nb.cols)
# # by quarter
# for(i in 1:length(plList)){
# plList[[i]] <- ggplot(datQ %>% filter(Prey == preys[i])) +
# geom_bar(aes(Year,eatenW,fill=Predator), stat="identity", position="fill") +
# scale_fill_manual(values=colPalette) +
# facet_wrap(~Quarter) +
# ggtitle(preys[i])
# }
# for(i in 1:length(plList)){
# postscript(paste(dirFigs,"/predators_on_",preys[i],".ps",sep=""))
# print(plList[[i]])
# dev.off()
# }
# by year
datY <- datQ20 %>%
group_by(Year,Predator,Prey) %>%
summarise(eatenW = sum(eatenW, na.rm=T))
#postscript(paste(dirFigs,"predators_on_eachprey.ps",sep="/"))
ggplot(datY) +
geom_bar(aes(Year,eatenW,fill=Predator), stat="identity", position="fill") +
scale_fill_manual(values=colPalette) +
facet_wrap(~Prey)+
xlim(1973,2020)+
theme(legend.position="bottom")
#dev.off()2017 SMS KeyRun
# ------------------------
# plot predators contributing to predation by prey
preys <- as.character(unique(datQ17$Prey))
plList <- vector("list",length(preys))
nb.cols <- 14
colPalette <- colorRampPalette(brewer.pal(12, "Paired"))(nb.cols)
# # by quarter
# for(i in 1:length(plList)){
# plList[[i]] <- ggplot(datQ %>% filter(Prey == preys[i])) +
# geom_bar(aes(Year,eatenW,fill=Predator), stat="identity", position="fill") +
# scale_fill_manual(values=colPalette) +
# facet_wrap(~Quarter) +
# ggtitle(preys[i])
# }
# for(i in 1:length(plList)){
# postscript(paste(dirFigs,"/predators_on_",preys[i],".ps",sep=""))
# print(plList[[i]])
# dev.off()
# }
# by year
datY <- datQ17 %>%
group_by(Year,Predator,Prey) %>%
summarise(eatenW = sum(eatenW, na.rm=T))
#postscript(paste(dirFigs,"predators_on_eachprey.ps",sep="/"))
ggplot(datY) +
geom_bar(aes(Year,eatenW,fill=Predator), stat="identity", position="fill") +
scale_fill_manual(values=colPalette) +
facet_wrap(~Prey)+
xlim(1973,2020)+
theme(legend.position="bottom")
#dev.off()